A non-equilibrium Monte Carlo approach to potential refinement in inverse problems 
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The inverse problem for a disordered system involves determining the interparticle interaction 
parameters consistent with a given set of experimental data. Recently, Rutledge has shown (Phys. 
Rev. E63, 021111 (2001)) that such problems can be generally expressed in terms of a grand 
canonical ensemble of polydisperse particles. Within this framework, one identifies a polydisperse 
attribute ('pseudo-species') a corresponding to some appropriate generalized coordinate of the sys- 
tem to hand. Associated with this attribute is a composition distribution p(a) measuring the number 
of particles of each species. Its form is controlled by a conjugate chemical potential distribution 
/i(o") which plays the role of the requisite interparticle interaction potential. Simulation approaches 
to the inverse problem involve determining the form of /i(o") for which p(o") matches the available 
experimental data. The difficulty in doing so is that /i(cr) is (in general) an unknown functional of 
p((t) and must therefore be found by iteration. At high particle densities and for high degrees of 
polydispersity, strong cross coupling between /i((j) and 'p{o') renders this process computationally 
problematic and laborious. Here we describe an efficient and robust non- equilibrium simulation 
scheme for finding the equilibrium form of ^[p(a)]. The utility of the method is demonstrated by 
calculating the chemical potential distribution conjugate to a specific log-normal distribution of 
particle sizes in a polydisperse fluid. 
PACS numbers: 61.20Ja, 02.70.Tt 



I. INTRODUCTION AND BACKGROUND 

Much of statistical mechanics is concerned with the 
task of deducing the macroscopic structure of a system, 
given a description of its microscopic interaction interac- 
tions. Sometimes, however, one is confronted with the 
reverse problem: that of determining the form of the 
interactions from knowledge of the structure. This so- 
called "inverse problem" jl| arises, for example, when 
one has scattering measurements of the structure factor 
of a molecular fluid and wishes to find the corresponding 
interparticle interaction potentials . Alternatively one 
might have obtained spectroscopic data for the distribu- 
tion of orientations of a certain bond in a molecular solid 
and wish to determine the corresponding bond orienta- 
tion potential 3]. 

One approach of longstanding for tackling inverse 
problems is Reverse Monte Carlo This simulation 
scheme seeks to minimize error estimators quantifying 
the difference between experimentally derived structural 
data (such as a radial distribution function) and the cor- 
responding simulation averages. The minimization pro- 
ceeds by replacing the Hamiltonian in the standard MC 
Metropolis update scheme with the error estimator func- 
tion, while the role of temperature is replaced by the es- 
timated degree of uncertainty in the experimental data. 
The method outputs a set of configurations, which are 
consistent with the experimental data and can be fur- 
ther analyzed, but provides no direct information on the 
underlying form of the interaction potential, the values 
of the thermodynamic fields, or fluctuation effects. 

Recently, Rutledge @ has suggested that ostensibly 
disparate inverse problems can be incorporated within 
a unifying theoretical framework by mapping them onto 
a generalized grand canonical polydisperse composition 
space. The key to achieving this is the identification of 



a generalized coordinate of the system, cr, on which the 
potential function of interest is defined. The problem 
may then be translated into the language of polydisper- 
sity by regarding cr as a continuous variable labeling the 
species of particles that comprise the system. In gen- 
eral, however, these polydisperse particles need not cor- 
respond directly to the real particles of the system, and 
are perhaps best regarded as 'pseudo particles' the pre- 
cise identity of which depends on the particular choice 
of cr and thus on the physical situation to hand. The 
reformulation is completed by assuming that the number 
density of pseudo particles of each a is free to fluctuate. 
Within the resulting open ensemble, the instantaneous 
count of particles of each species is measured by a com- 
position distribution p(cr), the ensemble-averaged form of 
which, p(cr), is directly controlled by the conjugate chem- 
ical potential distribution /i(cr) (see Appendix . The 
latter quantity plays the role of the potential function of 
interest. 

In practice a surprising variety of inverse problems can 
be cast in this form. For example, in a liquid crystal, 
one can identify cr with the orientation of a molecule (cf. 
0) and /i(cr) with an effective orientational potential. In 
an atomic crystal, cr can be identified with the phonon 
frequencies (5j and //(cr) to the Fourier transform of the 
interparticle potential. For a simple monatomic fluid, 
one can construct a "bond" picture, in which each of the 
N real particles is replaced by N non-interacting pseudo 
particles, all localized to the site of the real particle. Each 
pseudo particles is speciated according to its distance to 
one other nominated real particle and the associated ^i{o^) 
is simply the interparticle potential 0- 

In general, for systems described by two-body parti- 
cle interactions, a unique relationship is believed to exist 
between p(a) and /x(cr) 0,0, El- The simulation chal- 
lenge is then to determine the form of /t(cr) for which 73(cr) 
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matches some "target" form pj((7) corresponding to the 
available experimental structural data. The difficulty in 
achieving this is that /i(cr) is an unknown functional of 
p(cr), and thus must be determined iteratively via some 
refinement procedure. In the following section we sum- 
marize the state-of-the-art in this regard. 

II. ITERATIVE REFINEMENT SCHEMES AND 
CONVERGENCE ISSUES 

Refinement typically commences using some guessed 
form for the interaction potential /x(cr). This serves as 
input to a Monte Carlo simulation, in the course of which 
the form of p(cr) is obtained as a histogram. Once suffi- 
cient statistics for p((t) have been accumulated, the next 
and successive iterations modify /i(cr) to reduce the dis- 
crepancy between p(tT) and the target function pt{cr)- 
To achieve this Rutledge |5,] utilized an iteration scheme 
based on a zeroth order (i.e. non-interacting or ideal gas) 
approximation to the chemical potential distribution: 

^,^^<J)=\n{pia)) , (1) 

which can be viewed as the polydisperse generalization of 
the potential of mean force 11, 12]. Using this approxi- 
mation to initialize /^(cr), the iteration then proceeds ac- 
cording to 

M,:+i(a)=A..(a)-alnf^) , (2) 

where a < 1 is a step-size parameter, the value of 
which must be chosen sufficiently small to ensure that 
the method converges, but not so small that convergence 
is excessively slow. Rutledge applied this approach to 
calculate the effective interparticle interaction potential 
of a single component fluid consistent with a prescribed 
radial distribution function. For this particular purpose, 
his scheme is similar to the potential refinement method 
of Soper O. 

The efficiency with which inverse problems can be 
solved is clearly contingent on the convergence rate of 
the refinement procedure for /i(cr). Typically a strategy 
such as that of Eq.[21will be effective at low number den- 
sities of pseudo particles where the non-interacting ap- 
proximation retains some degree of accuracy. Difficulties 
arise, however, at high number densities and for high de- 
grees of polydispersity (ie. slowly varying forms of p(f )). 
In this regime, fJ.{(T) can differ greatly from p}^{a). More- 
over, strong cross coupling between /i(cr) and p(cr) implies 
that variations (within the refinement procedure) to the 
chemical potential at one value of a can significantly af- 
fect the entire form of p(it) . This can seriously impede 
convergence; the sole remedy being to utilize a small step 
size a, with a concomitant increase in the computational 
effort H. 



In view of these considerations, one might seek to 
adopt a more sophisticated refinement procedure in the 
expectation that it will promote faster convergence than 
Eq. 13 Such schemes have been developed and applied 
by Lyubartsev & Laaksonen flJl and Toth They 
employ large matrices of derivatives (obtainable via fiuc- 
tuation relations as ensemble averages) and designed to 
direct the minimization more reliably in the 'downhill' di- 
rection in the multidimensional parameter space of some 
cost function measuring the deviation of p(tT) from pticr). 
Whilst in our experience of gradient schemes (we have 
tried Powell's method, conjugate gradient methods and 
variable metric methods |lg), they do provide a more 
rapid approach towards the solution in the early iter- 
ations, we find that they do not converge as reliably 
(at least for high degrees of polydispersity) as the non- 
interacting approximation scheme of eq. [21 Instead we 
have found them to be prone to becoming trapped in lo- 
cal minima of the cost function. Although, empirically, 
this problem can be ameliorated by mixing iterations |l7| 
of the sophisticated scheme with simple ones of the form 
Eq.|21 such an ad-hoc approach is clearly neither theoret- 
ically well-founded nor aesthetically satisfying. 

Other recent work, carried out in the context of tar- 
geted polydispersity, has examined the utility of his- 
togram extrapolation in targeting a specific composition 
distribution pt (a) [llll9t. Histogram extrapolation |20j 
allows a histogram for p((t) accumulated at some /i(cr) 
to be reweighted to provide an estimate for p(tT) cor- 
responding to some other chemical potential distribution 
/i'(o'). To exploit the method, it can be embedded within 
a gradient-based scheme to minimize a cost function mea- 
suring the deviation of p(cr) and pticr) |19j . In practice, 
however, one finds that the extrapolation operates effec- 
tively only within a limited range of chemical potential 
space around the simulation state point. If one attempts 
to extrapolate too far from this point, spurious structure 
appears in the results and the method may not converge 
to the correct solution. This problem can be dealt with 
by constructing a series of intermediate targets that in- 
terpolate between the initial guess for p(cr) and the true 
target (see ref. for full details). 

In addition to the convergence issues already men- 
tioned, a further important matter is the statistical qual- 
ity of the measured estimate for 'p^{a). Unless the 'signal- 
to-noise' ratio associated with this measurement is suf- 
ficiently high prior to each iteration, the procedure will 
necessarily fail to converge. Often, however, one has no 
ready criteria for knowing when the statistical quality of 
the data is "good enough" . We note also in passing that 
related problems can arise when the statistical quality 
of the experimentally derived target distribution /9t(cr) is 
poor. Indeed, it has been reported that in such circum- 
stances the solution to which the refinement converges 
can depend on the initial guess for p,{a) or contain 
spurious features j22|. This would seem to suggest that 
'noise' in the target can engender additional local min- 
ima in the multidimensional parameter space on which 
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the cost function is defined. 

In this paper we describe a new iterative refinement 
scheme which is simple and straightforward to imple- 
ment, yet yields a robust and efficient procedure. Our 
method (which is inspired by, but is distinct from, the 
fiat energy histogram random walk algorithm of Wang 
and Landau '23]) is based on a non-equilibrium MC pro- 
cedure in which refinements to /i(fT) are implemented 'on- 
the-fly' during the simulation, not between iterations as 
in equilibrium schemes. Empirically we find that the 
method provides an estimate of the true solution right 
from the earliest iterations; subsequent iterations merely 
serve to reduce the statistical uncertainty in this esti- 
mate. This attractive feature promotes more reliable 
convergence than is found in more complicated gradient 
based schemes. In a comparative test, the new method 
is found to be significantly more efficient than the simple 
equilibrium iteration scheme described by eq. |2] 



Once C falls below some previously specified threshold 
value C*, A'(ct) can be saved and used as the starting point 
for the second iteration. This proceeds similarly to the 
first iteration, except that the value of the modification 
factor /3 is reduced, viz 

Pi+i=f3i/n, (5) 

with n some (small) integer |25l |. 

Successive further iterations steadily reduce f3 towards 
zero, thereby restoring detailed balance and driving /i(cr) 
towards its equilibrium limiting form. For sufficiently 
large i (i.e. sufficiently small Pi), ^ii{a) provides a good 
approximation to this limiting form. If desired, however, 
one can extrapolate fully to the equilibrium limit by per- 
forming a final run with /3 = 0, followed by application 
of histogram extrapolation to /i(cr) in order to 

match p((t) to pticr)- 



III. NON-EQUILIBRIUM REFINEMENT 
SCHEME 



We consider systems having a polydisperse attribute 
a varying continuously in the range < a < ac- The 
associated composition distribution p(cr) (defined in ap- 
pendix and its conjugate chemical potential distribu- 
tion are represented as histograms formed by binning the 
cr-domain into an integer number M of sub-intervals. A 
discussion of this binning strategy can be found in ref. 

ri9i. 

Given some prescribed target histogram pt(cr), the re- 
finement fJ.{a) proceeds iteratively using MC simulation. 
In general, no special bootstrapping measures are neces- 
sary for initializing fJ.{a) and, in the absence of a better 
guess, one can simply take p{a) = C \/ a, with C an 
arbitrary constant. During the course of an iteration, 
the chemical potential distribution is modified at regular 
short intervals (eg. every 10 MC sweeps). The modifi- 
cation continuously tunes fJ.{(T) such as to minimize the 
deviation of the instantaneous composition distribution 
p{a) from its target value: 



p{a) ~ ptja) 
Pt{cr) 



(3) 



where (3i is the modification factor for the ith iteration. 
Updating p{a) 'on-the-fly' in this manner quickly reduces 
differences between p{a) and its target |24l |. 

The criterion by which an iteration is deemed to have 
completed is based on the maximum value, across the 
range of cr, of the relative deviation of the average com- 
position distribution from the target form: 



max 



pti<y) 



(4) 



IV. RESULTS 

In order to gauge the efficacy of our scheme we have 
tested it on a challenging problem: namely that of de- 
termining the chemical potential distribution for a fiuid 
of polydisperse Lennard-Jones (LJ) particles having an 
extremely broad distribution of sizes. 

The interparticle potential between two particles la- 
beled i and j is given by 



with aij given by the additive mixing rule 



/ \ 12 













(6) 



(7) 



The target distribution of particle diameters a appearing 
in eq.Qwas assigned a log-normal form, described by the 
normalized shape function: 



[Ha /a) + (3/2) ln(l 



exp 



21n(l + 14^2) 



.(8) 



Here is the average particle diameter and W is the 
standard deviation in units of 'a. The associated tar- 
get composition distribution follows as /Ot((T) — Ptf{<^)i 
where pt is the target number density. 

For computational convenience, /(cr) was truncated at 
CTc = 12. This choice implies that the largest permit- 
ted particle size has a volume 1728 times that of the 
average particle size {a — 1). The simulations were per- 
formed within a periodic cubic box of side L = 90W and 
the temperature was set to T = 2.35 (in standard LJ 
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units ll^). The width and scale of the target polydis- 
persity distribution were set to = 2.5 and pt = 0.0152 
respectively, corresponding to a target volume fraction, 
77 « 13% where 

7? = rfffg . (9) 

Although this volume fraction would not be considered 
especially high for a pure fluid (the critical volume frac- 
tion for a monodisperse LJ fluid is about 16%), it is suf- 
ficient in the present context to ensure a high degree of 
cross coupling between p(cr) and /i((T). 

The grand canonical ensemble Monte Carlo (GCMC) 
algorithm employed has a Metropolis form [2^ and in- 
vokes four types of operation: particle displacements, 
particle insertions, particle deletions and particle resiz- 
ing; each is attempted with equal frequency. Specific to 
the polydisperse case is the resizing operation which en- 
tails attempting to change the diameter of a nominated 
particle by an amount drawn from a uniform random de- 
viate constrained to lie in some prescribed range. This 
range (maximum diameter step-size) is chosen to pro- 
vide a suitable balance between efficient sampling and 
a satisfactory acceptance rate at the prevailing number 
density. As regards the remaining types of moves, these 
proceed in a manner similar to the monodisperse case 
|2(tI |. except that for insertion attempts the new parti- 
cle diameter is drawn with uniform probability from the 
range a € [0, CTc]. 

Histograms of p(cr) and /i(cr) were formed by parti- 
tioning the range < a < ac into M — 120 equal inter- 
vals. The refinement procedure was initialized by setting 
/io(o') = V (T. The initial modification factor was set 
to Po = 0.01 (and halved at each successive iteration), 
while the threshold value of C at which an iteration is 
deemed complete (cf. Eq. 0J was set to (* — 0.4. With 
these parameters the simulation was found to converge to 
the correct chemical potential distribution in 7 iterations 
{/3j — 1.56 X 10^"*), requiring a total of approximately 
100 hours CPU time on an AMD 2000+ processor. The 
form of fi{(j) at the 1st, 3rd and 7th iterations are com- 
pared with the limiting equilibrium form in fig. |5J One 
sees that the principal effect of the reduction in /3 at each 
successive iteration is a decrease in the scatter of the so- 
lution; there seems to be no clear systematic dependence 
of the solution on the value of (3. This feature was con- 
firmed by independent runs at a number of other volume 
fractions and model parameters. 

Fig. El shows a snapshot of a typical configuration ob- 
tained using the equilibrium solution for fJ.{cr). The con- 
figuration contains some 1.1 x 10'' particles. Although the 
form of Pt (cr) (fig. ^ dictates that the largest particles 
occur with a probability of less than 0.5% that of the av- 
erage particle size, their far greater volume implies that 
they nevertheless occupy a sizeable fraction of the sys- 
tem volume. Herein lies the physical origin of the cross 
coupling between /)(cr) and p{a). Since the density (and 



chemical potential) of the small particles depends on the 
available volume, it must therefore also depend on the 
density (and chemical potential) of the larger particles, 
and vice-versa. 

Finally, for comparison, we have attempted to solve 
the same inverse problem using the equilibrium scheme 
of eq. 121 To initialize the chemical potential distribution, 
we took po{(y) = lii{pt{a)). Some preliminary tuning of 
the damping factor a and the run length of each itera- 
tion was found to be necessary in order for the method 
to converge. Trials with a — 1.0, a = 0.5 and a — 0.25 
failed to converge in a reasonable timescale. For a = 0.1 
the method did converge (in an oscillatory fashion) , tak- 
ing 93 iterations to reach the solution (compared to 7 
iterations for the non-equilibrium method), and in so do- 
ing consuming some 230 hours of CPU time. A sim- 
ilar run, starting with poicr) = V cr (as used in the 
test of the non-equilibrium approach) did not converge, 
thus demonstrating the sensitivity of the method to the 
quality of the initial guess. For this particular problem, 
the non-equilibrium scheme is therefore computationally 
more efficient by at least a factor of two. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have introduced a non-equilibrium it- 
erative potential refinement scheme for inverse problems. 
As a test, the method was used to determine the chem- 
ical potential distribution for a LJ fiuid having a wide 
distribution of particle sizes. We find it to be efficient, 
simple to implement and it requires no initial guess for 
the solution. 

The efficacy of the new method arises from its applica- 
tion of continuous small modifications to p{cr) through- 
out the refinement, rather than large-scale updates be- 
tween iterations as occurs in equilibrium schemes. The 
advantages of such an approach are apparent when one 
considers the high degree of cross coupling that can ex- 
ist between p{a) and p(o'), and the problem it implies: 
namely that a change to the chemical potential at any 
one value of a can substantially affect the whole compo- 
sition distribution. Accordingly it is desirable that (a) 
each individual modification to /^(cr) is kept small in or- 
der to moderate the magnitude of cross coupling effects, 
and (b) such deviations from the target distribution that 
do arise are immediately corrected by further modifica- 
tions to p{a). Our method achieves this by continually 
'tweaking' p{a) by small amounts to iron out differences 
between the instantaneous form of p{a) and the target. 

With regard to convergence properties, our findings 
indicate that even for the earliest iterations, the method 
yields (at least for current problem) an estimate for p{a) 
which, whilst statistically poor, does not appear to de- 
viate systematically from the true solution. As itera- 
tion proceeds, it "hugs" the solution ever more tightly, 
thereby allowing the degree of convergence to be judged 
solely on the basis of the statistical scatter of the solu- 
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tion. This remarkable feature would seem to promote 
reliable convergence to the true solution. One might ad- 
ditionally speculate that the intrinsically non-equilibrium 
character of the method generates "noise" that promotes 
escape from any local minima in the solution space, thus 
overcoming one of the apparent limitations of gradient- 
based refinement schemes. 

The main factor that we have found to influence the 
reliability and rate of convergence in our method is the 
choice of the initial value of the modification factor 
Po 113' Should this value be chosen to be too large 
(/3o > 0.05) then the method may require an excessive 
number of iterations to converge, or in extreme cases 
(as also occurs in equilibrium schemes ^J), not con- 
verge at all. If, on the other hand, it is chosen very 
small {Po < 0.005), then fewer iterations may be neces- 
sary overall, but the earlier iterations will take longer to 
complete because they require many more modifications 
steps. Within this window, the efficiency does not appear 
to vary by more than about a factor of 2, but there is cer- 
tainly scope for trial and error optimization. We found 
that this is achievable quickly (during the initial stages 
of the first iteration) simply by monitoring the time evo- 
lution of C (c.f . eq. . In cases where /3o is chosen too 
large for convergence, this quickly becomes evident dur- 
ing the first iteration because C does not decrease steadily 
towards C*- This permits early diagnosis of convergence 
problems, in contrast to equilibrium schemes, where a 



failure to converge is typically only apparent after a num- 
ber of iterations. 

Finally we note that although we have tested our 
method within the specific context of a polydisperse LJ 
fiuid, its potential applicability is more wide ranging. As 
described in sec.^ a diverse range of inverse problems can 
be mapped onto a grand canonical polydisperse composi- 
tion space 13 . It is simply a matter of identifying the ana- 
logue quantities to p{a) and //(cr); thereafter the method 
can be applied exactly as we have described. Thus, for 
instance, in a molecular liquid, a (partial) structure fac- 
tor g{r) might play the role of p{a) while the interatomic 
pair potential 0(r) plays the role of /i(cr). If applied to 
real experimental data in this way, it would be of inter- 
est to compare the robustness of our method with that of 
existing refinement schemes, particularly with regard to 
the sensitivity to experimental uncertainty in the target 
structural data and potential truncation effects p^ . 
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FIG. 1: The log-normal target distribution ptic) described 
in the text and displayed on a logarithmic scale. The distri- 
bution has standard deviation W = 2.5 and scale parameter 
pt = 0.0152. 
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FIG. 2; The convergence of the non-equilibrium refine proce- 
dure. Shown (bottom part) are the forms of lJ.i{o') at comple- 
tion of iteration numbers 1,3 and 7, together with the limiting 
equilibrium solution (solid line). Also shown for comparison 
is the potential of mean force f^^^{a) — ln(pt((j)). The up- 
per part of the figure shows (on the same scale) a difference 
plot measuring the deviation from the final solution at each 
iteration. 




FIG. 3: Snapshot configuration of A'^ = 11274 polydisperse LJ 
particles drawn from the log-normal size distribution shown 
in fig. El 



APPENDIX A: GRAND CANONICAL 
FORMULATION OF POLYDISPERSITY 

Here we provide a brief formal overview of the statis- 
tical mechanics of polydispersity. We consider a classi- 
cal fluid of polydisperse particles confined to a volume 
V = L"^. The system is assumed to be thermodynami- 
cally open, so that the particle-number distribution N(a) 
is a statistical quantity. The associated grand canonical 
partition function takes the form: 



N 



Zv = 



with 



N=0 ■ i=l 



dfi / dCTi ^ exp (-/JTYat ({r, o-})) 
V Jo ) 

(Al) 



N 



Hn ({f,a}) = $({f,a})-^M(aO 



(A2) 



Here N is the overall particle number, while (3 — 
{kBT)~^ and pl((j) are respectively the prescribed in- 
verse temperature and chemical potential distribution, 
{r, cr} denotes the configuration, i.e. the complete set 
(rl, (Ti), (r2, 0-2) • ■ • (rlv, ctn) of particle position vectors 
and polydisperse attributes. The corresponding configu- 
rational energy $ ({r, cr}) is assumed to reside in a sum 
of pairwise interactions 

N 

i<i=i 



where <j) is the pair potential. 
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The fluctuating particle number distribution is defined 

by 

N 

N{a) = Y,S{a-ai), (A4) 

i=l 

with a the continuous polydispersity attribute and N = 
J N{a)da. The associated composition distribution is 

p{a) = N{a)/V . (A5) 

The statistical behavior of p{a) is completely described 
by its probability distribution 



X exp(-/3W;v({f,a}))n5(p(a)-iV(a)/F) , 

(T 

from which the ensemble averaged form of pip) follows 
as 



-p{a) = \ p{u)vv\p{a)\dp{u) . (A7) 



